Optical Spectra of Solids Obtained by Time-Dependent Density-Functional Theory 
with the Jellium-with-Gap-Model Exchange-Correlation Kernel 



Paolo E. Trcvisanutto, 1, 2 Alcksandrs Terentjevs, 1 Lucian A. Constantin, 3 Valcrio Olcvano, 2,4 and Fabio Delia Sala 1,3 

1 National Nanotechnology Laboratory (NNL), Istituto di Nanoscienze-CNR, Via per Arnesano 16, 1-73100 Lecce, Italy 

2 European Theoretical Spectroscopy Facility (ETSF) 
3 Center for Biomolecular Nanotechnologies @UNILE, 
Istituto Italiano di Tecnologia, Via Barsanti, 1-73010 Arnesano, Italy 
4 Institut Neel, CNRS and UJF, Grenoble, France 
(Dated: October 29, 2012) 

Within the framework of ab initio Time Dependent-Density Functional Theory (TD-DFT), we 
propose a static approximation to the exchange-correlation kernel based on the jellium-with-gap 
model. This kernel accounts for electron-hole interactions and it is able to address both strongly 
bound excitons and weak excitonic effects. TD-DFT absorption spectra of several bulk materials 
(both semiconductor and insulators) are reproduced in very good agreement with the experiments 
and with a low computational cost. 



The theoretical description of the optical properties 
of materials by first principles calculations is one of the 
classical issues in solid state physics. After photon ab- 
sorption, the excitonic effects, driven by the electron 
and hole (e-h) interactions, are principal actors, render- 
ing the ab initio computational description demanding. 
The most successful approach is, so far, based on Many- 
Body Perturbation theory: the GW self energy accounts 
for electron-electron (e-e) many-body effects in the band 
structure calculations, whereas the Bethe-Salpeter Equa- 
tion is solved to introduce e-h interactions [l| . This accu- 
rate method has a limited applicability to large systems 
due to its high computational cost. 

Time-Dependent Density-Functional Theory (TD- 
DFT) is another theory for the exact treatment of ex- 
cited states. Similarly to ground state DFT, its main 
drawback is to properly approximate the unknown dy- 
namic exchange-correlation (xc) kernel f xc , which should 
account for both e-e and e-h interactions. From the very 
first attempt, the Homogeneous Electron Gas (HEG) 
model system had been employed in construction of 
xc kernels 0, 0|, and one of the most used approach 
was derived from local density approximation (LDA) 
xc potential v xc in the static regime, Adiabatic LDA 
(ALDA) 0. ALDA has been successfully applied to 
molecules and clusters, but it usually becomes inappro- 
priate in solids where the improvements with respect to 
the Random Phase Approximation (RPA) (with the triv- 
ial fxc = 0) are negligible. The introduction of short- 
range non-locality, like in a generalized-gradient approxi- 
mation (GGA) Q or in a non-LDA (NLDA) @ approach, 
does not improve with respect to LDA or RPA optical 
spectra. Main reason for this failure resides in the missed 
ultra-long range behaviour, 1 / |r — r'| or 1/q 2 in reciprocal 
space, absent in HEG and HEG-based kernels. On the 
other hand, the introduced empirical long range contri- 
bution (LRC) kernel, f xc = a/q 2 with a found re- 
lated to the screening, correctly reproduces e-h excitonic 
effects but it is limited to semiconductors and small gap 



insulators. Very recently, renewed interest in this issue 
has been boosted by two new approaches. In the first 
one, the Bootstrap kernel 0] extends the LRC kernel to 
a matricial form and proposes a heuristic form for the 
LRC weight in terms of the screening which has to be 
calculated self-consistently. This method has been suc- 
cessfully applied to a wide range of bulk materials. In the 
second work [l(|, starting from the meta-GGA xc func- 
tional they derive an f xc which presents an a/q 2 term, 
with a calculated from the meta-GGA parametrizations. 
Finally, the Nanoquanta kernel 0,111 adapted from the 
four-point Bcthe Salpctcr Equation kernel in the TD- 
DFT framework, achieves high accuracy, though without 
much less computational cost. 

In this letter, we propose a non-empirical static xc 
kernel based on the jellium-with-gap model (JGM) (ToT - 
17| . so relying on a well defined physics. In the spirit of 
DFT, this kernel is a density-functional f xc [n](E g ) and 
it also depends on the bandgap E g . The approach is ab 
initio in the same regards as the LRC, Bootstrap and 
Nanoquanta, once first principle calculations (like GW) 
are used to estimate good E g . We show that this kernel 
properly describes both weak excitonic effects in semi- 
conductors as well as bound excitons in ionic insulators, 
without requiring any frequency-dependence. It provides 
absorption spectra in good agreement with experiments. 
Moreover, this approach is as computationally expensive 
as standard ALDA. 

In linear- response TD-DFT, the central quantity is the 
density-density response function x(g,u;) (written in re- 
ciprocal space q and frequency to), which is calculated by 
the Dyson equation: 

X~ 1 (q,^ ! ) = Xo 1 (<1> 0J ) - f xc (q,oj) - v(q), (1) 

where xo{Qi^) is the independent-particle Kohn-Sham 
response function, and v(q) = 4ir/q 2 is the Coulomb po- 
tential. The dielectric function e is related to the density- 
density response via the relation e = 1 — v\- 
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For the HEG (or jcllium) model, which is the most 
important reference for simple metals, very accurate f xc 
kernels had been developed [l8l - [20j by fulfilling impor- 
tant exact constraints (e.g. compressibility and the third- 
frequency-moment sum rules) and using diffusion Monte 
Carlo input. 

Here we extend those parametrization to the case of 
the JGM model. In this Letter we restrict to the static 
uj = case and work in the adiabatic approximation. In 
the JGM model at density n and with energy gap E g , 
the RPA static dielectric constant (i.e. with f xc = 0) is 
El: 



JGM-RPA 
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Eq. ([2} is qualitatively different from the HEG coun- 
terpart, being finite for a non-zero band gap. Similar 
models are widely used to describe dielectric properties 
of semiconductors [2l[ . 

In analogy to Eqs. (20)-(21) of Ref. " 
for the JGM kernel (/i GM ) is approximated as: 



22 1 , a static kernel 



f^{q^0;n,E g ) 



-in 



1 



i 2 (4 



i) 



.3. 

nq 2 



(3) 

We then extend the kernel of Ref. [18| to the JGM im- 
posing condition ([3]) obtaining: 
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where k n , B{n) and C(n) are defined in Ref. [18j, being 
constructed from exact HEG constraints (k „ an d C(n)) 
and from HEG diffusion Monte Carlo data [111 (B(n)); 
kp is the Fermi wavevector. 

The JGM kernel has been constructed in order to sat- 
isfies the following properties: 

(i) f^ M (q; n, E g = 0)^ f™ G (q; »), where /™ G ( 9 , n) 
is the HEG static kernel defined by Eq. (12) of Ref. 
[l8l | , that is very accurate for the HEG correlation energy 
per particle at every wave vector q, and it recovers the 
accurate (exact) RPA long-range effective interaction; 

(ii) f^ M (q;n,E g ->• oo) -> —v(q), such that this ker- 
nel gives a vanishing correlation energy (E c ) in the limit 
of perfect insulators (E g — > oo). Moreover, due to the 
construction of B'(n,E g ) and C'(n,E g ), we observe by 
applying the so-called adiabatic-connection fluctuation- 
dissipation theorem [Lsl ]. that E c w constant/ E g + 
0(E~ 2 ), in the limit of large band gap. This is one of 



the most important exact constraints for the JGM which 
was derived from perturbation theory 171 ]; 
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Eq. ((7]) approaches to E 2 g jn (see Eq. j3])) for small E g 
(semiconductor case). Thus, in the q — > optical limit, 
we recover an LRC a/q 2 behaviour; 

(iv) lim^ 00 /i GM ( (Z ;n, J B s ) = -$-C'(n,E g ). This 
limit is unknown in the JGM, with the exception of the 
cases Eg — > and E g — > oo, for which the behaviour of 
our kernel becomes exact. We recall that both RPA and 
ALDA fail badly in this limit 0- 

In order to use the JGM kernel Eq. (fj| for real inho- 
mogencous bulk systems, we first consider its real space 
form / K J r GM (|r -r'\;n,E g ). We then use the DFT density 
n(r) in place of the constant density n, whereas as a first 
approximation (valid for bulk systems investigated in this 
work) we approximate E g to the fundamental band gap 
of materials; thus the kernel in the real space depends ex- 
plicitly on the position (i.e. /^ GM (|r — r'|, r; E g )). Then, 
we obtain the Fourier transform in the reciprocal space 
with the shape /i GM (|q + G|,G — G'\E g ), with q vec- 
tor in the first Brillouin zone, and G, G' reciprocal lat- 
tice vectors. The kernel still needs for symmctrization in 
G,G'. In previous works 0, @| the symmetrization was 
done in the real space. In analogy to the symmetrization 
adopted for the Coulombian term of the dielectric matrix 
231 ] , we perform the symmetrization in reciprocal space: 



/ KC (q, G, G'; E g ) = f:^ M ( ^|q + G||q + G'|, G-G'; E g ). 

(8) 

Eq. ((SJ) shows the following properties: 

(a) it reduces to Eq. ((3J in the homogeneous case. 

(b) The head, G=G'=0, element of the kernel matrix, 
in the optical limit approaches f xc (q — > 0,0, 0;E g ) w 
— (a)q~ 2 where 



a JGM (n(r),£ s )dr, 



(9) 



which differs from a JGM evaluated with the average den- 
sity fi. Thus, an explicit value for the LRC a parameter 
is provided as the mean value of a JGM in the unit cell. 
(a) depends on the bandgap and it is density-functional 
(i.e. it depends on the density inhomogeneity) . 

(c) In the q — > limit, the wings of the kernel matrix 
/«:(«, 0,G' ^ 0;E g ) = f xc (q,G ^ 0,0; E g ) cx l/|q + 
G||q+ GJ are 0(l/q) and the remaining elements are 
regular [3, [HI. 

(d) The ALDA-HEG terms are present in the kernel 
and they can play a role only at finite q, i.e. in energy- 
loss function calculations, where ALDA is already a good 
approximation [26]. 
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FIG. 1. (color online). Imaginary part of the macroscopic 
dielectric function (e2(o>)) for bulk Ge, Si and, GaAs. Red 
dashed line: RPA results with scissor operator shifting. Black 
solid line: TD-DFT with JGM kernel. Green dots are the 
experimental absorption spectra: Ge from ref. Si from 

ref. [33] , GaAs from ref . Same notation has been utilized 

in the next Figures. Ge and Si are obtained with relative 
Gaussian and Lorentzian broadening 0.02. Absolute Gaussian 
broadening 0.1 eV was used for GaAs. 
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FIG. 2. (color online). As in Fig. 1, but for SiC and Dia- 
mond bulk. Experimental absorption spectra: SiC from Ref. 
[3tJ and Diamond from Ref. [36fl . Absorption spectra are ob- 
tained with absolute Lorentzian broadening: SiC: 0.3 eV, and 
Diamond 0.2 eV. 



(e) The computational cost is very low, comparable to 
standard ALDA. 

The JGM kernel and its extension to inhomogeneous 
periodic systems (Eq. EJ), is based on several well defined 
and physical approximations which need to be tested in 
real situations. As a first test we applied Eq. ((8j) to 
calculate the optical absorption of several bulk systems. 

The absorption spectrum is calculated as the q — > 
limit of the imaginary part of the macroscopic di- 
electric function, $5e M (g, w), defined as e M (q,uj) = 
l/eQ^ G , =0 (c7, w). We have implemented the f xc kernel 
in the linear response TD-DFT dp code ^ . The DFT 
ground state calculations are performed utilizing ABINIT 
[28j . Like in all previous approaches 7 9], here we focus 
only on the e-h interaction, f^~ h , excitonic effects part 
of the full f xc , while we simulate the effect of the e-e, 
/|~ e , self-energy effects part by calculating xo with a 
scissor operator (so) bandgap corrected DFT-LDA elec- 
tronic structure, E so = E g — E^^, where E g is set to 
the experimental band gap [29| . The same parameter 
Eg is used into the JGM kernel. Instead of the exper- 
imental band gap, the GW one can be used as well, so 
that to refer to first principle calculations. Nevertheless, 
the problem to find the /|~ e kernel able to provide self- 
energy effects and the right E g (which is also an input 
of the JGM kernel), is once again postponed. Hereafter, 
the RPA results are meant to be scissor operator shifted. 

In Fig. [T] we show our calculated spectra for semi- 
conductors: Ge, Si and GaAs. These materials present 
weakly bound excitons nearby the onset of the contin- 
uum. Previous works Q have shown that e-h excitonic 
effects increase the oscillator strength of the lowest en- 
ergy part of spectra. 



Ge is the smallest band gap bulk system which we 
considered. Excitonic effects are modest, and the RPA 
calculation (dashed line) is already in good agreement 
with the experiment (green dots). Nevertheless, the TD- 
DFT JGM kernel (solid line) introduces improvements, 
enhancing the oscillator strength of the first peak at 2.2 
eV. 

On the other hand, in Si e-h interactions have impor- 
tant effects on the spectral weight of excitations and they 
determine a different absorption shape with respect to 
RPA. The JGM kernel increases the oscillator strength 
of the first feature at 3.4 eV where the RPA shoulder 
becomes a well defined peak, like in the experiment. 

Similarly to Si, in GaAs excitonic effects shift the spec- 
tral weight to lower energies. The JGM kernel moves the 
peak near 5 eV towards the experimental position and it 
increases the intensity of the structure at 3 eV therefore 
achieving a very good agreement with the experiment. 

More severe tests are Diamond and SiC (Fig. [2]) where 
the larger band gap E g increases the JGM kernel strength 
of excitonic effects. In SiC, the RPA peak at 8 eV is en- 
hanced and red shifted in the TD-DFT result, achieving 
an almost perfect overlap with the experimental peak. 
The higher energy peak at 9 eV is reproduced at the 
correct experimental position, though is quite overesti- 
mated. 

In Diamond, the position of the main RPA peak is cor- 
rectly red shifted by the TD-DFT, overlapping the exper- 
imental position, whereas the left shoulder is slightly too 
intensive. Overall, TD-DFT calculations show a smooth 
and realistic description of optical spectra of both SiC 
and Diamond. 

LiF and MgO are wide-gap insulators, with a strongly 
bound exciton occurring in the band gap 0, [l3| . They 
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TABLE I. The fundamental band gap E g [23], a/it from the 
best LRC fits, and (a) of Eq. (11). 
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FIG. 3. (color online). As in Fig. 1, but for LiF and MgO. 
Experimental absorption spectra: LiF from Ref. [3||, and 
MgO from Ref. [32j]. Here dotted-dashed magenta line is 
JGM TD-DFT only with diagonal f xc . Absorption spectra 
obtained with the absolute Gaussian broadening of 0.3 eV for 
both LiF and MgO. 



represent a severe workbench for any TD-DFT kernel. 
In Fig. |31 we show the calculated absorption spectra 
for LiF and MgO in comparison with the experiment. 
Here we report the results for both the full matrix JGM 
kernel (solid line), and also its G = G' diagonal part only 
(dot-dashed) , in order to show the effect of non-diagonal 
elements. The RPA calculations completely miss both 
the main peaks (LiF: 12.6 eV and MgO: 7.2 eV) and the 
resonant state features at higher energy. 

In LiF, the diagonal TD-DFT JGM kernel is able to 
conjure the main excitonic peak, albeit it is blue shifted 
by 0.5 eV with respect to the experiment and the inten- 
sity is overestimated. The second experimental peak at 
14.5 eV is well captured as shoulder, and it is joined in 
its final part. High energy features at about 20 eV are 



present, whereas the spurious BSE peak at 21 eV [13| 
here is absent. When inserting nondiagonal kernel ele- 
ments, we have further excitonic strength that red shifts 
the main peak up to 12.80 eV, improving the agreement 
with the experiment. Frequency dependence of f xc has 
been considered crucial for the correct description of the 
absorption spectra of wide band gap insulators where the 
strongly bound excitons are appearing 1^, 14 1. However, 
as shown in Ref. [§] and in this Letter, a static non-local 
kernel can also correctly address the e-h interactions in 
solid-state systems. 

In MgO, the strongly bound exciton is accompanied 
by high energy resonant peaks. JGM TD-DFT results 
display the strongly bound excitonic peak at 7.2 eV in 
agreement with the experimental findings. Likewise the 
LiF case, the height of the exciton is very intensive. The 
similarity with LiF concerns with the subsequent peak at 
10.3 eV: in JGM TD-DFT calculations, this becomes an 
artefact shoulder which joins the experimental results in 
the final part of the peak. The overall behaviour at high 
energies is in rather good agreement with experiments: 
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both energy and intensity peaks are well described. Non 
diagonal elements in the JGM kernel slightly decrease the 
intensity of the main peak at 7.2 eV, while the 13.50 eV 
peak is red shifted to 13.30 eV in better agreement with 
the 13 eV experimental result. 

Finally, we stress that the most important contribu- 
tion to the absorption spectra, is given by the head of 
the matrix kernel, that in the case of a static kernel, as 
the one presented in this work, is just an LRC-type a/q 2 , 
and all our results can be well reproduced by this simple 
(non-empirical) LRC-type kernel. In Table HI we show a 
comparison between JGM (a) values and previous a esti- 
mations: Results are well in agreement each other. This 
agreement begins to deviate for wider band gap bulk sys- 
tems where the inhomogeneity of the density enhances, 
resulting in much higher values of (a) for MgO and LiF. 

In conclusion, we have found that a simple, non- 
empirical and physical, static xc kernel derived from the 
jellium-with-gap model (JGM) provides accurate optical 
absorption spectra of a wide range of bulk solids (from 
metals to strong insulators). The physics of this kernel 
relies on the validity of JGM, in spite of HEG, to rep- 
resent real solids within TD-DFT. This kernel can be 
easily implemented in any solid-state code, and with a 
very reduced computational cost. In this work we tested 
bulk systems, but JGM kernel can be also applied to 
more complex systems, such as heterostructures or metal- 
organic interfaces, where the energy gap itself will depend 
on the position. The JGM kernel will correctly modulate 
the LRC contribution according to the local value of the 
gap. Further investigation will be devoted to energy-loss 
spectra at finite q. 
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